Variations

Here are some other equivalent ways to write this function.

The braces are not actually needed if the function only contains one statement.

fahr_to_kelvin <- function(temp) (temp-32) * (5/9) + 273.15

The body of the function within the {} can contain multiple lines of code, “statements”. Only the last line is returned.

fahr_to_kelvin <- function(temp) {
    kelvin <- (temp-32) * (5/9) + 273.15
    kelvin
}

return can be used to return from a function immediately. Further code will not be run. I recommend always using explicit return statements.

fahr_to_kelvin <- function(temp) {
    kelvin <- (temp-32) * (5/9) + 273.15
    return(kelvin)
    plot(1:10)
}

Statements in the function may be split over several lines, so long as it is clear that the statement is incomplete due to an unclosed bracket or an operator in need of a right hand argument.

fahr_to_kelvin <- function(temp) {
    kelvin <- 
        (temp-32) * (5/9) + 
        273.15
    return(kelvin)
}

Spacing and layout is largely for our own sanity.

   fahr_to_kelvin<- 
function(
                   temp){(temp 
  -32  )*(    5 
            /9 )+ 

273.15}

Beware! however of accidentally finishing a statement early.

fahr_to_kelvin_broken <- function(temp) {
    (temp-32) * (5/9) 
        + 273.15
}

fahr_to_kelvin_broken(212)

The first line in the body of the function is computed then discarded. The second line, + 273.15, is valid R so we don’t even get an error.

Finally, the computer doesn’t understand or care what names we give to functions, arguments, or variables.

charmander <- function(bulbasaur) {
    mew <- (bulbasaur-32) * (5/9) + 273.15
    mew
}

Composing Functions

Now that we’ve seen how to turn Fahrenheit into Kelvin, it’s easy to turn Kelvin into Celsius:

kelvin_to_celsius <- function(temp) {
    temp - 273.15
}

#absolute zero in Celsius
kelvin_to_celsius(0)

What about converting Fahrenheit to Celsius? We could write out the formula, but we don’t need to. Instead, we can compose the two functions we have already created:

fahr_to_celsius <- function(temp) {
    temp_k <- fahr_to_kelvin(temp)
    temp_c <- kelvin_to_celsius(temp_k)
    temp_c
}

# freezing point of water in Celsius
fahr_to_celsius(32.0)

This is our first taste of how larger programs are built: we define basic operations, then combine them in ever-larger chunks to get the effect we want. Real-life functions will usually be larger than the ones shown here–typically half a dozen to a few dozen lines–but they shouldn’t ever be much longer than that, or the next person who reads it won’t be able to understand what’s going on.

Function names and argument names are important cues to what code is doing.

Functions also help manage code complexity by cleaning up their contents after completion.

You might have noticed that all of these functions have an argument called temp. Why hasn’t this caused confusion and chaos? The answer is that arguments to a function and variables defined in the body of a function are local to a call to that function. Within the body of a function, R code can see variables from its own local environment, and variables from the global environment, but not variables local to other function calls. While running code R maintains a stack of local environments. When a call to a function returns, all of its local arguments and variables disappear.

Tip: There are several ways to see this in action.

• Scatter cats through your functions showing the values of arguments and variables at various points.

• Insert a call to browser() in one of your functions. This will pause execution and let you interact with the local environment in a function. Once in the browser, type “help” for a list of commands. Besides examining variables, you can step through the function from the point it paused with “n” or, at a finer grain, step into functions it calls with “s”. When you are done type “f”.

• You can debug an existing function using debugonce. This is like temporarily adding a call to browser() at the top of the function.

# debugging a function
debugonce(fahr_to_celsius)
fahr_to_celsius(212)

Running external software

Software outside of R is normally run from a “terminal” window, which is rather like the console in RStudio. You type in a command and the operating system runs it for you. The commands you type in are interpreted by something called a “shell”, because it acts like a shell around the operating system. For more information, see the Software Carpentry course on the Unix shell.

In RStudio, you can open a terminal using “Tools/Terminal/New Terminal”. Try this and type in:

uptime

The “uptime” command tells you how long the computer has been running for. R can give a command to the shell using the system function.

system("uptime")

Here is how we want to run FastQC for day 0. To R the command has no meaning beyond being a character string. It will be interpreted by the shell, which is external to R. The shell will then run FastQC for us.

system("fastqc --extract --outdir . r-progtidy-files/Day0.fastq")

In the “Files” tab of the bottom right pane in RStudio, you should see that some new files and a directory have been created. Click on “Day0_fastqc.html” to examine it. Also examine the file in the “Day0_fastqc” folder called “summary.txt”.

If you don’t have FastQC, or it failed to run for some reason: The expected output files can be found in the folder “r-progtidy-files/fastqc-output”.

Using a for-loop

Typing out the fastqc command for each day will get repetitive. Let’s use a for-loop to automate this task. The paste0 function lets us “paste” together character strings, so we can use that to construct the commands to run, like this:

# construct a command to run
day <- 0
command <- paste0("fastqc --extract --outdir . r-progtidy-files/Day", day, ".fastq")
command

We want to do this for each day, and then use system to run the resulting command:

# run fastqc on each of the files

days <- c(0,4,7,10,15,20)

for(day in days) {
    command <- paste0("fastqc --extract --outdir . r-progtidy-files/Day", day, ".fastq")
    cat("Running the command:", command, "\n")
    system(command)
}

Loading the summary.txt files

Note: If you weren’t able to run FastQC earlier, the expected output files can be found in the folder “r-progtidy-files/fastqc-output”. You will need to adjust the filenames appropriately in the R code below.

The summary.txt files are in “tab separated value” format. This is similar to comma separated value format we’ve seen in the introductory R day, but instead of using a comma it uses a special character called a tab to delimit columns. Tabs show up as variable amounts of space in a text editor. In R, they can be written "\t".

In base R we could use the function read.delim to read this file. It’s quite similar to read.csv. However, let’s use the Tidyverse package readr to load the file. To use readr we first need to load it with library.

library(readr)

We can now read the file like this:

read_tsv("r-progtidy-files/fastqc-output/Day0_fastqc/summary.txt")

Oh! There aren’t column headings in this file. We need to tell read_tsv this.

read_tsv("r-progtidy-files/fastqc-output/Day0_fastqc/summary.txt", col_names=FALSE)

We should tidy this up a bit before using it. Columns should have meaningful names, and the PASS/WARN/FAIL grading looks like it should be a factor:

filename <- "r-progtidy-files/fastqc-output/Day0_fastqc/summary.txt"
sumtab <- read_tsv(filename, col_names=FALSE)
colnames(sumtab) <- c("grade", "test", "file")
sumtab$grade <- factor(sumtab$grade, c("FAIL","WARN","PASS"))
sumtab

We expect to have to examine many FastQC reports in future. It will be convenient to have this as a function!

filename is going to be different each time we use this code, so it wants to be an argument to the function. The rest of the code goes in the body of the function. The returned value will be sumtab.

load_fastqc <- function(filename) {
    sumtab <- read_tsv(filename, col_names=FALSE)
    colnames(sumtab) <- c("grade", "test", "file")
    sumtab$grade <- factor(sumtab$grade, c("FAIL","WARN","PASS"))
    sumtab
}

load_fastqc("r-progtidy-files/fastqc-output/Day0_fastqc/summary.txt")

Applying the function

We now want to load each of the summary.txt files we created earlier. One perfectly valid way would be to use a for-loop. However here we’re going to use an “apply”-style function. “apply”-style functions are functions that call some other function repeatedly on subsets of a data set. R has a variety of functions with names such as apply and tapply, and several other functions where the idea is similar such as aggregate. apply itself is for use with matrices and produces a vector as output. Other functions in this family have different input and output types.

The flavour of apply we need here is lapply, which takes a vector and calls a function on each element in turn of that vector (i.e. the “subsets” here are simply individual elements). The result is a list containing the returned values, hence the name lapply.

First let’s work out the names of the files we want to load.

days <- c(0,4,7,10,15,20)
filenames <- file.path("r-progtidy-files/fastqc-output", paste0("Day", days, "_fastqc/summary.txt"))
filenames

Now we can use lapply to apply our function to each of the filenames.

sumtabs <- lapply(filenames, load_fastqc)

lapply is one of the base R functional programming tools. It stands for loop-apply and performs a specific class of loops without the user needing to write the loop structure. It can lead to clear code, when you are used to seeing it.

This example is equivalent to:

sumtabs <- list()
for (idx in 1:length(filenames)) {
  sumtabs[[idx]] <- load_fastqc(filenames[idx])
}

We now have a list of data frames. Lists are a little different from the vectors we usually work with. See section 20.5 in the “R for data science” book for a comprehensive explanation. Here are some ways to examine this object:

Basic differences are that lists can contain things of different types and different sizes. In this case, we have a list of data.frames, which are not a basic R type (a vector of something)

sumtabs
class(sumtabs)
length(sumtabs)
str(sumtabs)

When we print sumtabs on the console, it gives us a hint how to access individual elements: using the [[ ]] syntax:

sumtabs[[1]]
sumtabs[[2]]

This is an ungainly data structure, we would like to turn it into a single a single big data frame. The dplyr package provides a function to do this:

library(dplyr)
bigtab <- bind_rows(sumtabs)

bigtab
table(bigtab$test, bigtab$grade)

Base R has a similar idiom to do the same task. It is more flexible, but does less checking. I use it to bind vectors into a data.frame, which bind_rows refuses to do.

bigtab <- do.call(rbind, sumtabs)

Summary

What we’ve just done is:

1. Create a custom reader tool that sets up each file inside R

2. Use a looping tool to apply that function to a collection of file names

3. Transform the results produced by the looping tool into a data.frame that we can use for other purposes.

Improving load_fastqc

Our load_fastqc function will currently fail with an error if the file it is passed does not exist.

load_fastqc("nosuchfile.txt")

Depending on how we intend to use it, we might instead want the function to issue a warning and return NULL.

load_fastqc <- function(filename) {
    # Check arguments are sane
    if (!file.exists(filename)) {
        warning("No such file: ", filename)
        return(NULL)
    }
    
    # Load and tidy data
    sumtab <- read_tsv(filename, col_names=FALSE)
    colnames(sumtab) <- c("grade", "test", "file")
    sumtab$grade <- factor(sumtab$grade, c("FAIL","WARN","PASS"))
    sumtab
}

load_fastqc("nosuchfile.txt")

Checking input arguments before proceeding is a common pattern when writing R functions.